Asymptotic rate of quantum ergodicity in chaotic Euclidean billiards

The Quantum Unique Ergodicity (QUE) conjecture of Rudnick-Sarnak is that every eigenfunction phi_n of the Laplacian on a manifold with uniformly-hyperbolic geodesic flow becomes equidistributed in the semiclassical limit (eigenvalue E_n -> infinity), that is, `strong scars' are absent. We study numerically the rate of equidistribution for a uniformly-hyperbolic Sinai-type planar Euclidean billiard with Dirichlet boundary condition (the `drum problem') at unprecedented high E and statistical accuracy, via the matrix elements of a piecewise-constant test function A. By collecting 30000 diagonal elements (up to level n ~ 7*10^5) we find that their variance decays with eigenvalue as a power 0.48 +- 0.01, close to the estimate 1/2 of Feingold-Peres (FP). This contrasts the results of existing studies, which have been limited to E_n a factor 10^2 smaller. We find strong evidence for QUE in this system. We also compare off-diagonal variance, as a function of distance from the diagonal, against FP at the highest accuracy (0.7%) thus far in any chaotic system. We outline the efficient scaling method used to calculate eigenfunctions.Comment: 38 pages, 11 figures, version of Jan '06, in review, Comm. Pure Appl. Mat

Quasi-orthogonality on the boundary for Euclidean Laplace eigenfunctions

Consider the Laplacian in a bounded domain in R^d with general (mixed) homogeneous boundary conditions. We prove that its eigenfunctions are `quasi-orthogonal' on the boundary with respect to a certain norm. Boundary orthogonality is proved asymptotically within a narrow eigenvalue window of width o(E^{1/2}) centered about E, as E->infinity. For the special case of Dirichlet boundary conditions, the normal-derivative functions are quasi-orthogonal on the boundary with respect to the geometric weight function r.n. The result is independent of any quantum ergodicity assumptions and hence of the nature of the domain's geodesic flow; however if this is ergodic then heuristic semiclassical results suggest an improved asymptotic estimate. Boundary quasi-orthogonality is the key to a highly efficient `scaling method' for numerical solution of the Laplace eigenproblem at large eigenvalue. One of the main results of this paper is then to place this method on a more rigorous footing.Comment: 21 pages, 2 figures; preprint version of 7/4/0

Efficient numerical solution of acoustic scattering from doubly-periodic arrays of axisymmetric objects

We present a high-order accurate boundary-based solver for three-dimensional (3D) frequency-domain scattering from a doubly-periodic grating of smooth axisymmetric sound-hard or transmission obstacles. We build the one-obstacle solution operator using separation into P azimuthal modes via the FFT, the method of fundamental solutions (with N proxy points lying on a curve), and dense direct least-squares solves; the effort is O(N^3P) with a small constant. Periodizing then combines fast multipole summation of nearest neighbors with an auxiliary global Helmholtz basis expansion to represent the distant contributions, and enforcing quasi-periodicity and radiation conditions on the unit cell walls. Eliminating the auxiliary coefficients, and preconditioning with the one-obstacle solution operator, leaves a well-conditioned square linear system that is solved iteratively. The solution time per incident wave is then O(NP) at fixed frequency. Our scheme avoids singular quadratures, periodic Green's functions, and lattice sums, and its convergence rate is unaffected by resonances within obstacles. We include numerical examples such as scattering from a grating of period 13 {\lambda} x 13{\lambda} of highly-resonant sound-hard "cups" each needing NP = 64800 surface unknowns, to 10-digit accuracy, in half an hour on a desktop.Comment: 22 pages, 9 figures, submitted to Journal of Computational Physic

Unimodal clustering using isotonic regression: ISO-SPLIT

A limitation of many clustering algorithms is the requirement to tune adjustable parameters for each application or even for each dataset. Some techniques require an \emph{a priori} estimate of the number of clusters while density-based techniques usually require a scale parameter. Other parametric methods, such as mixture modeling, make assumptions about the underlying cluster distributions. Here we introduce a non-parametric clustering method that does not involve tunable parameters and only assumes that clusters are unimodal, in the sense that they have a single point of maximal density when projected onto any line, and that clusters are separated from one another by a separating hyperplane of relatively lower density. The technique uses a non-parametric variant of Hartigan's dip statistic using isotonic regression as the kernel operation repeated at every iteration. We compare the method against k-means++, DBSCAN, and Gaussian mixture methods and show in simulations that it performs better than these standard methods in many situations. The algorithm is suited for low-dimensional datasets with a large number of observations, and was motivated by the problem of "spike sorting" in neural electrical recordings. Source code is freely available

High-order discretization of a stable time-domain integral equation for 3D acoustic scattering

We develop a high-order, explicit method for acoustic scattering in three space dimensions based on a combined-field time-domain integral equation. The spatial discretization, of Nystr\"om type, uses Gaussian quadrature on panels combined with a special treatment of the weakly singular kernels arising in near-neighbor interactions. In time, a new class of convolution splines is used in a predictor-corrector algorithm. Experiments on a torus and a perturbed torus are used to explore the stability and accuracy of the proposed scheme. This involved around one thousand solver runs, at up to 8th order and up to around 20,000 spatial unknowns, demonstrating 5-9 digits of accuracy. In addition we show that parameters in the combined field formulation, chosen on the basis of analysis for the sphere and other convex scatterers, work well in these cases.Comment: 24 pages, 10 figure

How exponentially ill-conditioned are contiguous submatrices of the Fourier matrix?

We show that the condition number of any cyclically contiguous $p\times q$ submatrix of the $N\times N$ discrete Fourier transform (DFT) matrix is at least $$\exp \left( \frac{\pi}{2} \left[\min(p,q)- \frac{pq}{N}\right] \right)~,$$ up to algebraic prefactors. That is, fixing any shape parameters $(\alpha,\beta):=(p/N,q/N)\in(0,1)^2$, the growth is $e^{\rho N}$ as $N\to\infty$ with rate $\rho = \frac{\pi}{2}[\min(\alpha,\beta)- \alpha\beta]$. Such Vandermonde system matrices arise in many applications, such as Fourier continuation, super-resolution, and diffraction imaging. Our proof uses the Kaiser-Bessel transform pair (of which we give a self-contained proof), and estimates on sums over distorted sinc functions, to construct a localized trial vector whose DFT is also localized. We warm up with an elementary proof of the above but with half the rate, via a periodized Gaussian trial vector. Using low-rank approximation of the kernel $e^{ixt}$, we also prove another lower bound $(4/e\pi \alpha)^q$, up to algebraic prefactors, which is stronger than the above for small $\alpha, \beta$. When combined, the bounds are within a factor of two of the numerically-measured empirical asymptotic rate, uniformly over $(0,1)^2$, and they become sharp in certain regions. However, the results are not asymptotic: they apply to essentially all $N$, $p$, and $q$, and with all constants explicit.Comment: 24 pages, 4 figures. v3 slightly strengthens results via p not q in \sigma_1 bound, corrects minor typos (eg e\pi/4 inversions), updates abstract re K-B proo

Efficient high-order accurate Fresnel diffraction via areal quadrature and the nonuniform FFT

We present a fast algorithm for computing the diffracted field from arbitrary binary (sharp-edged) planar apertures and occulters in the scalar Fresnel approximation, for up to moderately high Fresnel numbers ($\lesssim 10^3$). It uses a high-order areal quadrature over the aperture, then exploits a single 2D nonuniform fast Fourier transform (NUFFT) to evaluate rapidly at target points (of order $10^7$ such points per second, independent of aperture complexity). It thus combines the high accuracy of edge integral methods with the high speed of Fourier methods. Its cost is ${\mathcal O}(n^2 \log n)$, where $n$ is the linear resolution required in source and target planes, to be compared with ${\mathcal O}(n^3)$ for edge integral methods. In tests with several aperture shapes, this translates to between 2 and 5 orders of magnitude acceleration. In starshade modeling for exoplanet astronomy, we find that it is roughly $10^4 \times$ faster than the state of the art in accurately computing the set of telescope pupil wavefronts. We provide a documented, tested MATLAB/Octave implementation. An appendix shows the mathematical equivalence of the boundary diffraction wave, angular integration, and line integral formulae, then analyzes a new non-singular reformulation that eliminates their common difficulties near the geometric shadow edge. This supplies a robust edge integral reference against which to validate the main proposal.Comment: 21 pages, 7 figures, revised version, to appear, J. Astron. Telesc. Instrum. Syst. (JATIS

High-order boundary integral equation solution of high frequency wave scattering from obstacles in an unbounded linearly stratified medium

We apply boundary integral equations for the first time to the two-dimensional scattering of time-harmonic waves from a smooth obstacle embedded in a continuously-graded unbounded medium. In the case we solve the square of the wavenumber (refractive index) varies linearly in one coordinate, i.e. $(\Delta + E + x_2)u(x_1,x_2) = 0$ where $E$ is a constant; this models quantum particles of fixed energy in a uniform gravitational field, and has broader applications to stratified media in acoustics, optics and seismology. We evaluate the fundamental solution efficiently with exponential accuracy via numerical saddle-point integration, using the truncated trapezoid rule with typically 100 nodes, with an effort that is independent of the frequency parameter $E$. By combining with high-order Nystrom quadrature, we are able to solve the scattering from obstacles 50 wavelengths across to 11 digits of accuracy in under a minute on a desktop or laptop.Comment: 22 pages, 9 figures, submitted to J. Comput. Phy

Validation of neural spike sorting algorithms without ground-truth information

We describe a suite of validation metrics that assess the credibility of a given automatic spike sorting algorithm applied to a given electrophysiological recording, when ground-truth is unavailable. By rerunning the spike sorter two or more times, the metrics measure stability under various perturbations consistent with variations in the data itself, making no assumptions about the noise model, nor about the internal workings of the sorting algorithm. Such stability is a prerequisite for reproducibility of results. We illustrate the metrics on standard sorting algorithms for both in vivo and ex vivo recordings. We believe that such metrics could reduce the significant human labor currently spent on validation, and should form an essential part of large-scale automated spike sorting and systematic benchmarking of algorithms.Comment: 22 pages, 7 figures; submitted to J. Neurosci. Met

Explicit unconditionally stable methods for the heat equation via potential theory

We study the stability properties of explicit marching schemes for second-kind Volterra integral equations that arise when solving boundary value problems for the heat equation by means of potential theory. It is well known that explicit finite difference or finite element schemes for the heat equation are stable only if the time step $\Delta t$ is of the order $O(\Delta x^2)$, where $\Delta x$ is the finest spatial grid spacing. In contrast, for the Dirichlet and Neumann problems on the unit ball in all dimensions $d\ge 1$, we show that the simplest Volterra marching scheme, i.e., the forward Euler scheme, is unconditionally stable. Our proof is based on an explicit spectral radius bound of the marching matrix, leading to an estimate that an $L^2$-norm of the solution to the integral equation is bounded by $c_dT^{d/2}$ times the norm of the right hand side. For the Robin problem on the half space in any dimension, with constant Robin (heat transfer) coefficient $\kappa$, we exhibit a constant $C$ such that the forward Euler scheme is stable if $\Delta t < C/\kappa^2$, independent of any spatial discretization. This relies on new lower bounds on the spectrum of real symmetric Toeplitz matrices defined by convex sequences. Finally, we show that the forward Euler scheme is unconditionally stable for the Dirichlet problem on any smooth convex domain in any dimension, in $L^\infty$-norm.Comment: 37 page